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ABSTRACT 

The structure and stability of flames in dusty mixtures is 
investigated. The presence of the dust leads to significant transport of 
energy by radiation and the fundamental goal of the analysis is to explore to 
what extent this displaces the classical non-hydrodynamical stability 
boundaries of the plane deflagration. An approximate description of the 
radiative transport permits analysis for arbitrary values of both the Planck 
length and the Boltzman number. It is shown that the pulsating /traveling- 
wave instability usually associated with values of Lewis number (Le) bigger 
than 1 is strongly enhanced by the presence of radiation and can be present 
even if Le < 1. On the other hand radiation tends to suppress the cellular 
instability normally associated with values of Le less than 1. The latter is 
consistent with preliminary experimental observations of Abbud-Madrid 
and Ronney. 


’This research was supported by the National Aeronautics and Space Administration under NAS1-19480 
while the authors were in residence at the Institute for Computer Applications in Science and Engineering 
(ICASE), NASA Langley Research Center, Hampton, VA 23681. 
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Introduction 

We are concerned with plane flames supported by mixtures 
that contain a large number of small inert particles (dusty gases). 
The particles substantially influence the radiative transport and 
make emission and absorption energetically significant. Partial 
motivation for the study comes from recent experimental 
observations which suggest that in combustion fields of this kind 
cellular instability (the Turing instability associated with values of 
Lewis number less than 1, [1]) can be suppressed [2]. A discussion of 
stability and the effect of particle loading on the location of the 
stability boundaries in the wave-number - Lewis-number plane is 
the central contribution of this paper. 

The presence of particles and radiative transport presents 
serious technical obstacles to an analytical treatment, and these must 
be overcome by judicious modeling. We start by assuming that the 
mass loading is small and the only impact of the particles is on the 
radiative transport. The bulk particle density is( 4/3) 7id 3 Np (N is 
the number density, d the particle diameter, p the substantive 
particle density) and the radiation effect is controlled by Nd^ so that 
our treatment is formally valid in the limit N —> », d — » 0,Nd 2 fixed. 
(Then Nd 3 0) 

A strategy for dealing with radiative transport which has been 
usefully employed by Joulin and his colleagues in a number of 
important studies (e.g. [3],' [4]) is to assume that the Boltzman 

number is small, but we wish to avoid the restriction here since it is 
of limited applicability. Instead we start with the differential - 
equation approximation (the Eddington approximation) 
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L v( L V-q) = 3q+4* L VT 4 . (1) 

Since the nonlinearities here present serious difficulties we 
approximate the Planck length L by a constant and linearize the 
emission term, replacing Eqn. (1) by 

L 2 v(v- q) = 3q + 1 8 Tc L V T (2) 

where T c is a constant characteristic (emitting) temperature. Then q 
is irrotational so that we may introduce a generating function \|/ so 
that 

q = V\|/ . ( 3 ) 

Equation (2), after a single integration, is then equivalent to 

L 2 V 2 V - + 18Tc L(T -T f ) (4) 


where Tf is the remote cold-mixture temperature. 

The other ingredients of our model are familiar and of proven 
merit. They are embodied in the equations 

9 -E/2FT* 

p(Y t + uY x ) = pD V t - Be 8(n) 

(5) 

9 9 -E/2FT. 

pCp(T t + uT x ) = X V T -V \|t + QB e 8(n) 

Here Y is the mixture mass fraction, T the temperature, and u is the 
steady flame speed so that the unperturbed flame is fixed. Reaction 
is modeled by a delta-function whose strength is an Arrhenius 
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function of the flame-temperature T*; n is the distance measured 
along the flame-sheet normal that points towards the burnt gas. And 
for the purposes of the stability analysis we adopt the constant- 
density model so that hydrodynamic effects are discarded. 

Boundary conditions are: 

x — > -oo T» Tf , Y-> Yf , — > 0 (6) 

We wish to construct the steady planar solution of the system (4) - 
(6) and to examine the linear stability of this solution to arbitrary 
perturbations. To this end it is convenient to introduce a flame-fixed 

coordinate system by replacing x by s where 

s^x-f(y,) ( 7 ) 

and the flame-sheet is located at s = 0. Equations (5) for s ^ 0 are 
then 


p [Y t + ( u -F t ) Y J = pDV 2 Y , 
p C p [T, +(u -F t ) TJ =XV 2 T -V 2 y 


( 8 ) 


where 


v 2 he (i +f 2 ) a 2 /as 2 + a 2 /ay 2 - z= y a 2 /ayas -F yy a /as (9) 


The jump conditions at the flame-sheet are: 


M - 0 ,[vj = 0 
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[T] - 0 ,[Y] - 0 


( 10 ) 


(l +F/j[Tj = -(QBA.) exp( - ERT.) , (l + F y =) [Yj -(B<>D) exp( - ERT.) . 


Steady Solution 

The steady solution is constructed by seeking solutions for 
which 

F = 0 , T =T S (s) , Y =Y s (s) , \| f =\|/®(s) ( 1 1 ) 

Consider the cubic 

a(coL) 3 -(coL) 2 -(3a + p)((ol_) +3 = 0 (12) 


where a =ULMC P , p =ieT c 4 /MC p T c . 


(13) 


a is the familiar length X/M C p divided by the Planck length; p is a 
Boltzman number. 

Equation (12) only has real roots, two of which are positive 
(coi , 0)2), one negative (0)3). In addition we define 0)4 b y 

£04 =pu /pD = M /pD . (14) 

Then, without using the connection conditions (10), the steady 
solution has the form: 


CD., S C0 2 s 


s < 0 T s =T f + C 1 e 1 + C 2 


e 


V s = Ci (M Cp -Xa>i)e 1 -C2(M C p -A,a> 2 )e 


C 0 o s 


(15) 
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CD. S 

Y s - Y f + C 4 e 

0)„ S 

Sj>J2 T s =Ta + C 3 e 

V s = C 3 (M C p -X^e^ 5 (16) 

Y = 0 . 

Here T a = Tf +CYf/C p is the adiabatic flame temperature, the 
temperature at s— fixed by global energy conservation. 

The various constants {Ci} are determined by use of the 

connection conditions (10 a-d) whence 

C 4 = Yf (17) 

and, with the definitions 

Zj = aL coi , Di =Ci(T a -Tf)- (18) 

we have 

Di + D 2 - D 3 = 1 

Di zi + D 2 Z 2 - D 3 z 3 = 1 (19) 

Di z 1 2 + D 2 Z 2 2 - D 3 z 3 2 = 1 


Equations (19) have solution 
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D^ = (z 2 - )(Z3 - )(z 2 -Zi) '\z 3 -Zi) - 1 
D 2 = (z-) - J(Z 3 - )(Zi -z 2 ) "\z 3 -Z 2 ) - 1 


D 3 - -(zi - )(z 2 - >(zi -z 3 ) -1(zg -Z 3 ) - 1 (20) 

In order to determine the mass flux M (i.e. the flame speed) it is 
convenient to define M c by 

M 0 - M(p - () (21) 

and since T. = T a at this limit, Eqn. (lOe) is 

[T«| - (MoCp/X) (Ta-Tf) exp[E/2FT a - E / 2R] . (22) 

It follows that 

M /M 0 = exp((E/2RTa)D 3 [D 3 + T a /(T a -T,)] i . (23) 


D 3 is defined by the roots Zi , z 2 and 2q which, in turn, are defined by 
a and (3, both of which depend on M. Thus Eqn. (23) is an implicit 
formula for the flame-speed. 

The limit a -* C 

It is of interest to examine the realistic limit a C when the 
Planck length is much larger than the diffusive scale (i.e. the 
Reynolds number based on L is large). Then the temperature 
distribution is 
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£jsJ2 T s = Tf + (T a -T,)[exp((MCp A.) 4 + fjp . 4 ] . 

SL2.Q T s = T a + (T a -Tf) f .(p , 4 (24) 


where 


f±(p .4 = p(p 2 + li '’'expj [(-1 p ± 1 Vp 2 + 12)/ l] sj 


On the scale s = O(L) there is a radiative preheat and post he at zone 
separated by a thin region (s = O (\/MC p )) in which the temperature 
increases from 


to 


T, = T, +(T a -Tf) p(p 2 + l4 ■’ (25) 

T. * T a +(T«-T,) p(p 2 + l4 '' ,2 (26) 


The flame-speed is the adiabatic flame speed corresponding to 
a cold mixture of temperature Ti so that 


M /M 0 



(EfiRTaJptp+Tada-TfC'lp 2 



(27) 


In order to discuss this formula it is convenient to introduce 
the virtual Boltzman number 

p 0 = 16 Tc 4 /M 0 C p T c 


(28) 
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Using CFU/air adiabatic flame calculations of Giovangigli and Smooke 
[6] and with T c equal to the adiabatic flame temperature we can plot 
variations of (3 0 with equivalence ratio, and these are shown in Fig. 1. 
This defines a realistic range of values of p 0 • 

Figure 2 shows variations of M/M 0 and P with p 0 determined 
from Eqns. (27), (28). The enhancement in the burning rate was 
described in ref. [3] in the case of small p (p = 0(RT a /^) and we call it 
the Joulin effect. It is analogous in its physical origins to the 
enhancement described by Weinberg in so-called excess enthalpy 
flames [5]. Here the temperature at the reaction zone is raised above 
the adiabatic flame temperature by radiative preheating, whereas in 
Weinberg's flames the increment is generated by providing a heat 
conduction path superior to that of the gas. 

The limit a - » 0 , /J — > ° © , a/J = 0( 1 

The formulas (24) and (27) are not uniformly valid in p. 
Indeed, there is a distinguished limit a -» 0 ,p -4 ,aP = 0(1) In 
this limit 


coiL ~ (2tx) _1 [l +Vl + 4cP ] , C 02 L ~ & /aP 

(29) 

C 03 L ~ (2a) - 1 [l Vi + 4tp] 

Thus the radiative preheat zone has thickness 0(l_ 2 MC p A.) 
whereas the postheat zone has thickness 0(A. /MC p ). The 
temperature rises slowly from Tf to T a ; increases rapidly from T a to 
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T. ~T a + (T a -Tf)(l + 4<p) • 1 ' 2 (30) 

and then decreases rapidly to T a . The burning rate is given by the 
formula 

M/M 0 = exp{ (E/2RT a )[l +T a (T a -Tf) ‘Vi + ] _1 } . (31) 

It is well known that in the limit of infinite activation energy 
the reaction zone is unstable if the temperature falls off towards the 
burnt gas on the same scale (O (X/MC p )) as it falls off towards the 
fresh gas [7]. This is why, for example, the so-called premixed-flame 
domain of diffusion flame burning, for which the burning rate is a 
decreasing function of Damkohle number (a portion of the middle 
branch of the familiar S-shaped response), is unstable. Thus the 
steady solution described in this section may not be physically 
realizable.* 

a Not Small 

When a is not small we return to Eqn. (23) in order to calculate 
the burning rate. Some curves are shown in Fig. 2. It is noteworthy 
that for high particle loadings, when a is not small, and for large 
values of the virtual Boltzman number, there is little augmentation of 
the burning velocity. This result should be compared with the 


* 


The reaction zone instability occurs on a fast-time scale with the combustion 
field beyond the reaction zone unchanged. It is not captured by the stability 
calculations of the present paper. 
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saturation observed experimentally [2]. Note that in an experiment 
Po is varied by changing the equivalence ratio, and a 0 can be varied 
with Po fixed by changing L (i.e. by changing the particle loading). 

Stability Analysis 

We seek infinitesimal perturbations of the steady solution by 
writing 


T = T s + eTo (s)e 


iky +a>t 


Y . Y s + eY 0 (s)e 


iky +cot 


iky +tot 

¥ “¥^ + £¥ 0 ( s )® 

_ ^ iky +(ot 

F - e6e 
where e— > 0. 

The coefficient functions satisfy the equations: 


pC p [co T 0 + u T 0 ' -8coT s ] -X [t 0 ”-k 2 T c + 8k 2 T s ] 

- [¥ 0 "-k 2 ¥ 0 + 8k 2 \K s ] 

p [coY o+uY o -8co Y s ] - pD[Y 0 "-k 2 Y 0 + 8k 2 Y s ] 


L 2 [y 0 -kVo + Sk 2 ^] = 3v]/ o + 1 &T c 3 LT c 


Solutions are: 

T, s Yj> s to.s to, s 

s < 0 T 0 = Ai e +A26 +8G1 e + 8G26 


(32) 


( 33 ) 
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Y 0 =Ce^ + 8 G 7 e t04S (34) 

Y 1 s s co 1 s co Q s 

y 0 = B-j e + B 2 © + 6G4© + SGs© 

s > 0 T 0 -A :i e : ‘ S +A i e ,e + 8 G 3 e’* S 

Y 0 = 0 , (35) 

Y3 S Y4 s ^ CO- S 

V 0 * B3 © + B4 © + 5 G 6 © 

Here y , y 2 »Y 3 an ^Y 4 are the roots of the quartic 

( 7 L ) 4 -a ^(yL ) 3 -(yLf[2{kLf + 3 +a -1 (p +coLu *)] + (yL)cc ^[(kL ) 2 + ^ 

+ (kL) 4 +(kL) 2 (3 +a _1 (p +coLu ■)) + 3tx VoLu _1 = 0 (36) 

where Re ^ ,y 2 )>C and Re ^.Y^O The constants Gi through G7 
are known (particular solutions) and C, Ai, A 2 , A 3 , A 4 are unknown. 
The constants {Bi } are related to { Ai } by 

Bi - 18T c 3 L [L 2 ( Yj 2 -k 2 ) . (37) 

Moreover (i is given by the formula 

P .1 M(pD)-' + l[M 2 (pD)' 2 + 4(k 2 + pco/pD)] v2 . (38) 


Thus there are six unknown constants (A i , A2 , A3 , A4 , C ,5) . 
The jump conditions (10), when linearized, are: 
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[T 0 ] - 0 ,[ Vo ] - 0 , [vj - 0 , [Y 0 ] - 0 , 


(39) 

[t„] = - BTa -T() To /2Rlf aL , [y „] == Elf , M T 0 /2Rlf pD . 


These provide six relations amongst the six constants and so define a 
linear homogeneous system which only has a non-trivial solution for 
certain values of the eigenvalue co. In this way, stability boundaries 
in the wave-number, Lewis-number plane can be constructed. 

Representative stability results are shown in Fig. 3 for an 
activation energy 0 of 15. Plotted are neutral stability boundaries 
for (a) a = 1, (b) a = 0.5, (c) a = 0.1, and (d) a= .02, in the Lewis 
number (l_e ='k /pD C p ) - scaled wave-number (ak) plane, and for 
various values of the Boltzman number (P). 

The right stability boundary (denoted by dashed curves in Fig. 
3) corresponds to a pulsating or traveling wave instability (Im co * 
0), one that is not accessible for common mixtures when p = 0. But 
radiative effects push this boundary sharply to the left, to commonly 
achieved values of Le. Figure 4a shows the locus of the intersection 

of this boundary with the line ak = c in the Le-P plane for c = 
.001, 0.1, 0.5, and 1.0, when a = 1.0. The shift to the left (smaller Le) 
with increasing p is eventually reversed. This is no different from 
the effect of burner-induced heat losses on the stability boundary in 
which displacement to the left is followed in due course by reversal 
of the entire curve [8]. 

The right stability boundary is not shown in Fig. 3c (a = .02). 
For small values of a the roots are extremely sensitive to variations 
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in P and there are serious convergence difficulties. Sensitivity to a 
parameter which is only roughly defined (because of the linearized 
transport equation) discouraged us from examining this situation 
more carefully. If it is to be done it should be done in the context of 
a more accurate model than the one considered here. 

The analogy with burner-induced losses carries over to the left 
stability boundary (denoted by solid curves in Fig. 3) which 
corresponds to the familiar cellular flames and for which co = 0. 
There is a strong tendency to suppress the instability by 
displacement of the curve away from Le = 1. Ronney [2] has 
reported that lean methane/air flames are smoothed by the addition 
of dust, in agreement with these results. Cellular flames in dustless 
mixtures are usually unsteady, a consequence of the fact that the 
maximum Lewis number on the stability boundary is then achieved 
when k = 0 (see, for example, the numerical solutions of the 
Kuramoto-Sivashinsky equation in [9]). This characteristic does not 
in general survive the addition of dust. For example, when a =0.1, P 
= 5 the right most point on the boundary is located at a k 21 0.42. 

Under this circumstance we anticipate a steady cellular structure on 
a scale defined by this special value of the wave number. The same 
effect is seen when burner-induced heat losses are accounted for [8], 
as seen in stationary polyhedral flames [10], [11]. 

Concluding Remarks 

We have shown that the introduction of significant radiative 
transport by the addition of dust to combustible mixtures can have a 
significant effect on the thermal-diffusive stability boundaries. The 
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pulsating instability can be strongly enhanced and the cellular 
instability suppressed. Suppression of cellular instabilities in lean 
methane/air mixtures has been reported by Ronney [2]. 

The radiation model that we have adopted is a crude one, 
although it appears to retain the essential physics. One of its 
characteristics is that it leads naturally to the definition of a 
Boltzman number that contains the numerical factor 16 rather than 4 
(cf. Eqn. (13b)). As a consequence too much significance should not 
be attributed to the precise values of (3 that appear in the discussion 
and in the figures. The order-of-magnitude is correct however. 
Should more detailed experimental data become available in the 
future it might be appropriate to carry out a more accurate analysis, 
together with more detailed exploration of the stability boundary 
locations. 
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Fig. 1 Virtual Boltzman number vs. Equivalence Ratio for lean 
methane/air flames, calculated using the flame temperatures 
and speeds obtained numerically by Giovangigli and Smooke [6]. 
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RIGHT BRANCH 
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Fig. 4 Locus of the point Re(co) = 0, ak = c in the Lewis number - p plane. 
4a. Right branch, 0 = 15 , a = 1.0, T a /Tf = 5. 

4b. Left branch, 0 = 15 , a = 1.0, T a /Tf = 5. 
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